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Abstract 

A  new  algorithm  for  computing  quantile  regression  estimates  for  problems  in  which  the 
response  function  is  nonlinear  in  parameters  is  described.  The  nonlinear  /  \  estimation 
problem  is  a  special  (median)  case.  The  algorithm  is  closely  related  to  recent 
developments  on  interior  point  methods  for  solving  linear  programs.  Performance  of  the 
algorithm  on  a  variety  of  test  problems  including  the  censored  linear  quantile  regression 
problem  of  Powell  (1986)  is  reported. 


Keywords:  Quantile  Regression,  Nonlinear  Regression,  Linear  Programming,  Interior 

Point  Algorithms,  Nonlinear  Programming 


* 
Department  of  Economics,  University  of  Illinois,  Champaign,  Illinois,  61820.    This  work  was  partially 
supported  by  NSF  grant  SES  91-21776. 


1.      Introduction 

About  a  century  ago  Edgeworth  observed  that  methods  of  estimation  based  upon  minimiz- 
ing sums  of  absolute  residuals  could  be  far  superior  to  least-squares  methods  under  non- 
Gaussian  error  conditions.  Laplace  had  drawn  similar  conclusions  a  century  earlier.  See  Stigler 
(1986)  on  this  early  history.  But  computation  of  l[  -estimators,  even  for  linear  regression, 
remained  a  major  impediment  to  applications  until  the  emergence  of  the  simplex  algorithm  for 
linear  programming  in  the  1940's.  Papers  by  Chames,  Cooper  and  Ferguson  (1955)  and  Wagner 
(1959)  provided  a  foundation  for  modern  algorithms  for  linear  l\  -regression  by  Barrodale  and 
Roberts  (1973),  Bartels  and  Conn  (1980)  and  others.  These  algorithms  are  readily  extended  to 
linear  quantile  regression,  Koenker  and  Bassett  (1978),  of  which  /  \  -regression  is  an  important 
(median)  special  case. 

The  current  state  of  algorithms  for  nonlinear  quantile  regression  is  far  less  satisfactory. 
Certainly  nothing  comparable  to  the  venerable  Gauss-Newton  algorithm  for  nonlinear  least 
squares  problems  has  emerged.  Despite  a  flurry  of  interest  by  prominent  numerical  analysts  in 
the  1970's  and  early  1980's,  see,  e.g.,  Osborne  and  Watson  (1971),  Murray  and  Overton  (1981) 
and  Bartels  and  Conn  (1982),  occasional  applications  of  nonlinear  quantile  regression  have 
relied  on  the  Nelder  and  Mead  (1965)  algorithm  and  other  generic  optimization  methods.  An 
excellent  statement  of  the  current  state-of-the-art  is  provided  in  the  thesis  of  Busovaca  (1985). 

In  contrast,  the  statistical  theory  of  nonlinear  quantile  regression  has  developed  rapidly  in 
recent  years.  Powell  (1986)  has  emphasized  its  value  in  the  analysis  of  censored  and  truncated 
responses.  Asymptotic  theory  for  independent  errors  has  been  developed  by  Oberhofer(1982), 
Dupacova(1987),  and  Powell(1991).  Theoretical  developments  by  Weiss  (1991)  and  White 
(1991)  have  stressed  applications  to  time-series  analysis.    Applications  of  Horowitz  and  Neu- 


mann  (1987),  Chamberlain  (1990),  and  others  have  demonstrated  its  value  in  applied 
econometrics. 

In  this  paper  we  will  describe  a  new  approach  to  the  computation  of  nonlinear  quantile 
regression  estimators  based  on  recent  interior  point  methods  for  solving  linear  programs.  In  the 
next  section  we  review  interior  point  methods  for  strictly  linear  problems.  Section  3  describes 
our  approach  to  nonlinear  problems,  and  Section  4  describes  our  computational  experience. 

2.      Interior  Point  Methods  for  Linear  Programs 

In  this  section  we  provide  a  brief  discussion  of  interior  point  methods  for  solving  strictly 
linear  programs  including  the  linear  quantile  regression  problem.  Our  exposition  will  follow 
closely  that  of  Vanderbei,  Meketon,  and  Freedman  (1986)  and  Meketon  (1986).  We  should 
emphasize  that  in  our  experience  interior  point  algorithms  for  linear  quantile  regression  do  not 
appear  to  be  competitive  in  efficiency  terms  with  existing  simplex  method  algorithms.  See 
Koenker  and  d'Orey(1987)  for  a  description  of  a  simplex  based  algorithm  for  linear  quantile 
regression.  However,  unlike  simplex  based  methods  they  do  appear  to  offer  a  natural  extension 
to  nonlinear  problems.  Thus  a  clear  understanding  of  the  linear  case  is  an  essential  first  step  in 
our  exposition. 

2.1.   A  Canonical  LP 

Consider  the  equality  constrained  linear  program 

min  {c'co|coe  Q={coe  R£,  Au  =  b  }  }  (2.1) 

where  R+  denotes  the  positive  orthant  of  Rn.  Given  a  feasible  point  in  the  interior  of  the  con- 
straint set,  co  e  int  (H),  interior  point  methods  proceed  in  two  steps.  First  we  transform  coordi- 
nates to  reposition  co  so  it  is  centered  relative  to  the  set  Q.  Then  a  (projected)  gradient  step  is 


taken  toward  the  boundary  of  Q.  Repeating  this  process  brings  us  arbitrarily  close  to  a  solution, 
and  a  stopping  criterion  is  eventually  invoked. 

To  flesh  out  this  brief  description,  let  D  =  diag  (co)  and  consider  the  transformation 

co  -»D-1co 

We  have  D_1co  =  1„,  an  ^-vector  of  ones,  so  the  transformation  D  has  the  effect  of  centering  co 
relative  to  the  orthant  boundaries  of  CI.  Correspondingly,  we  may  define  A  =  AD  and  c  =  Dc.  In 
the  transformed  coordinates  we  wish  to  move  in  the  gradient  direction  -c,  but  to  preserve  feasi- 
bility we  should  instead  project  c  onto  the  null  space  of  A  to  insure  that  the  equality  constraints 
are  satisfied. 

Let  c  denote  this  projection,  i.e., 

c={I  -A'{A'ATlA)c 

Clearly,  c  is  a  direction  of  descent;  and  we  now  move  toward  the  boundary  of  Q  in  this  direction. 
Let 

.A 

ot  =   max    {e{c  } 

;  =  1 n 

where  et  is  the  i     unit  basis  vector  for  R".  For  some  fixed  r\  e  (0,  1),  consider 

co  <-  co  -  (r\/a)Dc 
which  defines  a  sequence  of  iterations  co^  =  r(co^).  Since  at  each  iteration 

c'co*+1  =c,(ak-(r\/a)c'Dc  =c'coj.-(r|/a)||c||, 
we  expect  to  see  an  improvement  in  the  objective  function  at  each  iteration. 

Proposition.    If  c  <  0  the  problem  (2.1)  is  unbounded,  unless  c  =  0  in  which  case  every  co  e  Q  is 


optimal.  Otherwise,  the  problem  is  bounded  and  the  sequence  {c'co*  }  is  stricdy  decreasing. 

Proof.  Since  the  proof  of  this  proposition,  found  in  Vanderbei,  Meketon  and  Freedman  (1986),  is 
both  elementary  and  revealing  we  repeat  it  here  for  the  sake  of  completeness.  If  c  =  0,  there 
exists  a  vector  z  such  that  c  =  A'z,  hence  Dc  -  DAz  and  since  co  €  int  (Q)  it  follows  that  c  =  Az. 
But  then  for  any  co  e  Q., 

c'o)  =  z'Aa)  =  z'b 

which  is  independent  of  co,  establishing  that  c'co  is  constant  on  all  of  Q.  Next  consider  c  <  0. 
Note  that 

c'co!  =  c'co  -  7(co)c'c  =  c'co  -  7(co)  ||  c  II 2  (2-4) 


where  y(co)  =  r\/a.  Since  c  <  0, 


cop  =  co  -  pDc 


is  feasible  for  any  p  >  0  and 


c'cOp=c'co-p  ||  c  \\i 


implies  that  c'cop  ->  -oo  as  p  ->  <».  Finally,  if  c  >  0,  then  since  7(co)  >  0  and  c'co!  <  c'co  follows 
from  (2.4),  establishing  that  the  step  is  a  direction  of  descent.         ■ 

2.2.  Linear  /  \  -regression 

In  the  linear  model 

yt  =Xj'$  +  Uj  i  =  1,   •  •  •  ,  fl, 

as  noted  in  the  introduction,  the  /  j  -estimator  of  (3  which  minimizes 


R(b)=  £|y,  -Xi'b  |. 

1=1 

may  be  formulated  as  a  linear  program.  The  dual  problem  may  be  written  as 

max  [y'd  \d  e  Q  =  {  d  e  [-1,  If ,  X'd  =  0  }  } 

where  y  is  the  ^-vector  of  responses,  and  X  is  the  n  x  p  design  matrix.  To  solve  the  dual  problem 
we  proceed  as  before,  except  that  the  centering  is  slightly  altered  to  accommodate  the  altered 
form  of  Q..  For  any  initial  feasible  point  d,  e.g.,  d  =  0,  following  Meketon  (1986),  set 

D  =  diag  (min  { \+dt,  \-d,  }  ). 

In  the  transformed  coordinates  D~ld  the  projected  gradient  is 

Du  =  (I  -DX\X'D2XylX'D)Dy  =  D(y  -Xb) 

where  b  =  (X'D~X)~[D~y.  Note  that  as  in  the  former  case  the  transformation  has  the  effect  of 
centering  the  point  d  in  the  feasible  set  Q..  Now  let 

ex  D  u      -et  D  u 

a  =  max  { max  { ,  }  } 

\+di         l-di 

and  again  for  r\  e  (0,  1)  we  take  the  step 

d  <r-  d  +  (r\/a)D2u. 

Note  the  change  in  sign  since  we  are  now  maximizing.  The  iteration  sequence  d^+i  =T(d/c)  in 
the  dual  vector  implicitly  defines  a  corresponding  primal  sequence  with 

bk  =  (X'DiXY'X'Dly. 

As  Meketon  notes,  the  duality  theory  yields  a  natural  stopping  criterion.  Since 

y'dk  ^Z  \yi-*i%\ 


with  optimality  if  and  only  if  equality  holds,  it  is  reasonable  to  stop  iterating  when  the  differ- 
ence between  the  dual  and  primal  values  is  less  than  a  specified  tolerance. 

2.3    Linear  Quantile  Regression 

If  we  replace  the  (symmetric)  / 1  -criterion  with  an  asymmetric  linear  criterion  so  we  minim- 
ize 

*e(*)=£ Pete -*.■*) 
1=1 

pe(w)  =  (6  -  I(u  <  0))u,  we  obtain  the  regression  quantiles  of  Koenker  and  Bassett  (1978).  The 
dual  problem  is  now, 

max  [y'd\d€  Q={de  [0-1,  0]",  X'd  =  0]  ) 

This  leads  to  an  algorithm  identical  to  the  / 1  special  case  except  that  now 

D  =diag(min(Q-d;,  1-0+4)) 

and 

'  e{Dc        -e{Dc 


a  =  max  (max 


0-^'   1-0  +  4 


3.      Nonlinear  Quantile  Regression 

To  extend  these  ideas  to  the  case  of  nonlinear  response  functions  we  begin  by  considering 
the  nonlinear  l{  problem 


min  £|)5(0I  (3.1) 


where,  for  example, 


fi(0  =  yi-fo(xi,0. 

As  noted  by  El  Attar,  et  al  (1979)  a  necessary  condition  for  t  to  solve  (3.1)  is  that  there  exists  a 
vector  d  e  [-1,  1]"  such  that 

J(t*)'d  =  0  (3.2) 

A**)'*  =  Z  !//('*)  I  (3.3) 

where/  (r)  =  (/•(*))  and/(r)  =  0/(0%) 

Thus,  as  proposed  by  Osborne  and  Watson  (1971),  one  approach  to  solving  (3.1)  is  to  solve 
a  succession  of  linearized  l\  problems  minimizing 

Il/i(0-i,(0'5|  =ll/--/5|h, 

choosing  a  step  length,  X,  at  each  iteration,  by  line  search  in  the  resulting  directions  5.  The 
difficulty,  as  we  see  it,  with  this  approach  is  not  only  that  we  must  expend  the  effort  to  solve  an 
1 1  linear  program  at  each  iteration,  but,  perhaps  more  significantly,  the  resulting  directions  may 
actually  be  inferior  to  directions  determined  by  incomplete  solutions  to  the  sequence  of  linear- 
ized problems. 

Let  t  be  the  value  of  the  parameter  at  the  current  iteration,  and  consider  the  dual  problem 

max  [f'd  e  [-1,  If,  J'd  =  0}.  (3.4) 

If  the  model  were  linear  so 

f{s)=f{t)-K{s-t) 

for  some  fixed  matrix  K,  then  a  solution  can  be  found  by  applying  Meketon's  algorithm  to  find 
d    to  solve  (3.4),  computing 


5*  =  (K'D2K)-lK'D2f. 

where  D  =  diag  (min  { 1-d; ,  1+d,  })  and  setting  t  =  t  +  5  .  When  /  is  nonlinear  there  is  no 
longer  a  compelling  argument  for  fully  solving  (3.4)  at  each  iteration,  indeed,  in  our  experience 
only  a  few  iterations  to  refine  the  dual  vector  is  preferable.  In  the  version  of  the  algorithm  we 
have  implemented  to  conduct  the  tests  reported  in  the  next  section  we  take  two  dual  steps 
between  successive  updates  of /and  J.  A  detailed  description  of  the  algorithm  is  now  provided. 

3.1.  Dual  Step 

For  any  feasible  d  in  the  interior  of  the  constraint  set  of  (3.4)  we  refine  d,  following  Meke- 
ton,  as  follows.  Let 

D  =diag  (min  (1-d,-,  \+d,  }) 

s  =D2(I  -J(J'D2JTlJD2)f. 

d  <r-d  +  (r\/a)s 

where 

a  =  max  {  max  {$,7(1  -  dfi,  -5,7(1  +  d,-)  }  } 

and  r\  e  (0,  1)  is  the  constant  chosen  to  insure  feasibility.  Following  Meketon,  we  use  r\  =  .97. 
Updating  D,  s,  and  the  new  d  continues  the  iteration.  This  process  is  embedded  in  a  sequence  of 
primal  iterations  in  which  we  update /and  J  as  follows. 

3.2.  Primal  Step 

The  dual  step  yields  the  primal  direction 


b  =  (J'D2jylJ'D2f' 

which  we  explore  by  line  search.  Our  current  implementation  uses  Brent's  (1973)  algorithm 
from  the  PORT3  library  Fox  (1984).  Updating  we  have 

t  <-  t  +  \*h 

where  X   =  argmin  \\f(t  +  A.S)||i,  and  we  then  update  /and  J.   However  before  returning  to  the 

dual  step  we  must  adjust  the  current  d  so  that  it  is  feasible  for  the  new  value  of/.  This  is  accom- 
plished, somewhat  naively,  by  projecting  the  current  d  onto  the  null  space  of  the  new  /,  i.e. 
d  =  (/  -J(J'J)~lJ')d  and  then  shrinking  it  to  insure  that  it  lies  in  [-1,  1]",  so 

d  <-  d/(max  { |  di  \  }  +  e) 

for  some  tolerance  parameter  e  >  0.  Obviously,  when  the  problem  is  linear,  so  J  is  fixed,  this 
"adjustment"  is  nugatory  since  d  is  already  in  the  null  space  of  7,  and  the  algorithm  is  essentially 
like  Meketon's. 

3.3.  Stopping 

The  algorithm  currently  terminates  when  the  new  iterate  fails  to  improve  the  objective 
function  by  a  specified  tolerance.  Exploration  of  alternative  stopping  rules  is  a  topic  for  future 
research. 

3.4    Related  Literature 

Gill,  Murray,  Saunders,  Tomlin,  and  Wright  (1986)  and  Bayer  and  Lagarias  (1991)  have 
recently  pointed  out  the  close  relationship  of  "projected  Newton  barrier"  methods  (see  Fiacco 
and  McCormick  (1965))  and  interior  point  methods.  Algorithms  closely  related  to  the  one 
described  above  could  presumably  be  formulated  employing  logarithmic  barrier  functions  in  the 
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dual  vector  d. 

3.5.  Quantile  Regression 

As  in  the  case  of  the  linear  problem  the  generalization  of  the  /  j  problem  to  other  quantiles 
is  straightforward  involving  only  a  modification  of  the  constraint  set  [-1,  1]"  to  [0-1,  Q]n  for 
some  9  e  (0,  1). 

4.      Numerical  Experience 

In  this  section  we  describe  our  computational  experience  with  a  variety  of  test  problems. 
To  facilitate  comparison  with  existing  results  in  the  literature  we  have  chosen  problems  from 
Busovaca(1985)  and  Wormersley(1986).  We  focus  exclusively  on  the  I  \  case  since  there  are  no 
comparable  results  in  the  literature  for  other  quantiles.  The  problems  used  are  described  in 
detail  in  Appendix  A.  We  have  attempted  to  investigate  all  of  the  problems  reported  on  by 
Busovaca,  however  in  several  cases  we  were  unable  find  a  complete  description  of  the  problem. 
The  problem  taken  from  Wormersley  is  included  to  explore  the  important  special  case  of  piece- 
wise  linear  reponse  functions  which  arise  in  Powell's(1986)  formulation  of  the  quantile  regres- 
sion problem  for  censored  data. 

In  Appendix  B  we  provide  explicit  versions  of  our  algorithm  in  both  S,  Becker,  Chambers 
and  Wilks(1988),  and  Gauss.  All  of  the  reported  tests  were  carried  out  in  S  on  a  Sun  3/50.  Note 
that  the  line  search  algorithm  in  the  S  and  Gauss  versions  are  different.  To  implement  a  simple 
version  of  the  Osborne  and  Watson(1971)  algorithm  in  S  we  employ  the  S  function  11  fit 
which  does  l\  regression  using  the  Barrodale  and  Roberts(1973)  algorithm.  The  S  function 
lsf  it  carries  out  the  corresponding  weighted  least  squares  computations  for  the  interior  point 
algorithm. 
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A  summary  of  our  experience  on  the  test  problems  appears  in  Tables  5.1  and  5.2..  For 
Wormersley's(1986)  censored  regression  problem  (Problem  1)  our  version  of  the  interior  point 
algorithm  converges  to  Wormersley's  reported  solution.  However,  it  should  be  noted  that  the 
solution  to  this  problem  is  notoriously  nonunique.  Busovaca's  algorithm  cannot  be  employed  on 
Problem  1  due  to  the  fact  that  the  Hessian  of  the  response  function  is  identically  zero  almost 
everywhere.  The  remaining  problems  are  all  taken  from  Busovaca,  and  generally  our  interior 
point  solutions  correspond  closely  to  his.  In  Problems  7  and  13  there  are  small  discrepancies 
favoring  Busovaca;  in  Problem  9  there  is  a  larger  descrepancy  favoring  the  interior  point 
method.  Results  for  our  implementation  of  the  Osborne  and  Watson  algorithm  are  somewhat 
less  satisfactory.  It  fails  completely  on  Problems  11  and  12,  performs  poorly  in  Problems  1,  4b, 
and  13,  but  does  slightly  better  than  the  interior  point  method  on  Problem  5.  All  three  algorithms 
fail  for  Problem  4a  which  is  highly  degenerate  at  the  specified  initial  point.  At  an  alternative 
starting  point,  the  interior  point  algorithm  performs  well. 

5.      Some  Concluding  Remarks 

We  have  described  a  simple  approach  to  computing  quantile  regression  estimates  for  prob- 
lems with  nonlinear  response  functions.  The  approach  is  based  on  recent  developments  on  inte- 
rior point  methods  for  linear  programming,  but  may  be  viewed  as  a  variant  of  well-known  itera- 
tively  reweighted  least  squares.  While  the  algorithm  seems  to  perform  well  on  a  variety  of  test 
problems,  there  is  considerable  room  for  improvement.  Handling  rank  deficiency  in  the  model 
Jacobian  is  critical.  Alternative  stopping  criteria  also  should  be  explored. 
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Figure  5.1 
Algorithmic  Performance  on  Several  Test  Problems 


Interior  Point 

Algorithm 

Osbome-Watson  Algorithm 

Busovaca  Algorithm 

Example 

Starting 

Optimal 

Number  of 

Optimal 

Number  of 

Optimal 

Number  of 

Point 

Objective 

Iterations 

Objective 

Iterations 

Objective 

Iterations 

l.Wormersley 

(0,0) 

3.032544 

3 

5.234825 

3 

NA 

2.Bard 

(1.1,1) 

0.1243555 

6 

0.1243383 

5 

0.1243406 

13 

3.Beale 

(1,0.1) 

0.2928905e-07 

6 

0.154761  le-06 

6 

0.3695488e-05 

8 

4.a.Biggs 

(1,1,1,1,1,1) 

F 

F 

F 

b.Biggs 

(1,8,2,2,2,2) 

0.0 

11 

0.7289559 

45 

NA 

5.Brown&Dennis 

(25,5,-5,-1) 

903.3648 

29 

903.2406 

113 

903.2343 

2 

6.E1- Attar  5.1 

(1,2) 

0.47042 

10 

0.4704267 

6 

0.4704247 

8 

7.E1-Attar5.2 

(1.1.1) 

7.902733 

13 

7.904731 

22 

7.894227 

5 

S.Madsen 

(3,1) 

1.0 

16 

1.000010 

11 

1.000002 

13 

9.0sborne  1 

(0.5,1.5,-1,0.01.0.02) 

0.0293912 

14 

0.0293914 

10 

0.8203727 

55 

1  O.Osborne  2 

(1.3,0.65,0.65,0.7,0.6, 
3,5,7,2,4.5,5.5) 

F 

F 

F 

11. Powell 

(3,-1,0,1) 

0.1272e-07 

15 

F 

0.29039e-08 

3 

12.Rosenbrock 

(-1.2,1) 

0.0 

15 

F 

0.506642e-06 

51 

B.Watson 

(1,1.1,1) 

0.6487749 

25 

1.278432 

5 

0.6018584 

24 

H.Wood 

(0,0,0,0) 

0.0 

7 

0.000003 

10 

0.0 

25 

|  See  Appendix  A  for  a  detailed  description  of  the  test  problems. 
F  indicates  the  algorithm  failed  to  meet  convergence  criteria  for  the  problem. 
NA  indicates  results  are  not  available  for  this  entry. 


Table  5.2 
Optimal  Points  for  Several  Test  Problems 


Example 

Interior  Point  Algorithm 

Osborne-Watson  Algorithm 

Busovaca  Algorithm 

l.Wormersley 

-6.74166.4.59299 

-7.29804,4.74178 

NA 

2.  Bard 

0.10095.  1.52545.  1.97182 

0.10094,  1.52516,  1.97211 

0.10094,  1.52513,  1.97214 

3.Beale 

3.0.5 

3,0.5 

2.99999,  0.49999 

4a.Biggs 

F 

F 

F 

b.Biggs 

1,  10,  1,5,4,3 

1.82143,81.94978.2.27882 

NA 

5.Brown&Dennis 

-9.70273.  11.74096,  -0.442O4 

-10.0227.  11.91354,-0.44026 

-10.2236,  11.90843,  -0.45804 

0.55827 

0.55823 

0.58032 

6.E1-Attar5.1 

2.84250,  1.92018 

2.84250.  1.92018 

2.84250,  1.92018 

7.E1- Attar  5.2 

0.53558,  -0.00139,  0.02871 

0.53148,-0.00004,0.02751 

0.53606,0.0,0.00319 

8.Madsen 

0.0,0.00016 

0.0,  0.0022 

0.0,  -0.00205 

9.0sbome  1 

0.37706,2.19244,-1.72549 

0.37706.2.19246,-1.72552 

1.06716,  1.80257,-1.80731 

0.01332,0.02129 

0.01332,0.02129 

0.00345.  0.00109 

lO.Osbome  2 

F 

F 

F 

11. Powell 

0.7268e-O4.  -0.7268e-05 

F 

0.5588e-08,  -0.3725e-09 

0.1162e-04.0.1163e-04 

0.1250e-08,0.1716e-08 

12.Rosenbrock 

1.0,  1.0 

F 

0.99999,  0.99999 

13.  Watson 

-0.37526,  1.14089 

-0.23584,  1.03241 

-0.44271,  1.19321 

-0.42239,  0.39683 

-0.22747,0.41384 

-0.47676,  0.38449 

14.Wood 

1.0,  1.0,  1.0,  1.0 

1.0,  1.0,  1.0,  1.0 

1.0,  1.0,  1.0,  1.0 

Appendix  A 

Test  problem  1  (Wormersley,  1986) 

The  results  of  temperature  accelerated  life  tests  on  electrical  insulation  in  40  mo- 
torettes  are  recorded  in  Table  A.l.  This  data  is  originally  from  Schmee  and  Hahn  (1979). 
Ten  motorettes  were  tested  at  each  of  four  temperatures.  Testing  was  terminated  at  dif- 
ferent times  at  each  temperature.  The  model  used  to  fit  the  data  is 


1000.t2 

loqiQii  =  xi  +  — h  e, 

y  (T  + 273.2) 


where  H  is  the  failure  time  and  T  is  the  temperature. 


Table  A.l  :  Data  for  motorettes  example 


150 


Test  temperature  T°C 
170  190 


200 


Failure  times  H  in  hours 


1764 

408 

408 

2772* 

408 

408 

3444 

1344 

504 

3542 

1344 

504 

3780 

1440 

504 

4860 

5196 

8064 

5448 

1680 

528 

10  units 

3  units 

5  units 

5  units 

Termination  time  H 


*  Wormersley  gives  the  second  failure  time  at  170°  as  2722,  but  his  results  are  consistent 
with  the  values  recorded  here  from  Schmee  and  Hahn. 


At  each  temperature  there  is  an  upper  bound  H  (the  time  at  which  testing  was 
stopped)  on  the  observed  failure  times,  so  the  algorithms  of  the  observed  failure  times  are 
given  by 

.     /,        r-r                 1000x2 
min     loq in H ,  x,  H — — (-  e 

V  (T  + 273.2) 

Test  problem  2  (Bard,  1970) 

(  Ui 

Ji\x)  =Vi-  [xi  + 


V{X2  +  WiX3 

where  i  =  1,2,  •  •  •  ,  15,  Ui  =  i,  V{  =  16  —  i,  ivt  =  mm(w,-,Vj),  and 


i 

yt 

i 

Vi 

i 

Vi 

1 

0.14 

6 

0.32 

11 

0.73 

2 

0.18 

7 

0.35 

12 

0.96 

3 

0.22 

8 

0.39 

13 

1.34 

4 

0.25 

9 

0.37 

14 

2.10 

5 

0.29 

10 

0.58 

15 

4.39 

Test  problem  3  (Beale,  1958) 

fi(x)  =  Vi  -3l(l  -*2)» 

where  z  =  1,2,3,  y1  =  1.5,  y2  =  2.25  and  y3  =  2.625. 
Test  problem  4  (Biggs,  1971) 


fi(x)  =  x3exp(-ttxi)  -  x4exp(-t1x2)  +  x6exp(-tlx5 )  -  y2, 
where  i  =  1,  •  •  • ,  13,  £z  =  (O.l)z  and 

yi  =  exp(-ti)  —  bexp(  —  10tl)  +  3exp(  —  4t{). 

Test  problem  5  (Brown  and  Dennis,  1971) 

fi(x)  =  (xi  +  tix2  -  exp(ti))2  +  (x3  +  i4«in(<,-)  -  cos(i2))2, 

Test  problem  6  (El- Attar  5.1,  1979) 

/i(ar)  =  *}+z2-10 

/2(x)  =  xi  +  x\  -  7 
f3(x)  =  x\-x\-l 

Test  problem  7  (El- Attar  5.2) 

fi(x)  =x\  '+  x\  +x\  -  1 
/2(x)=x;  +  xl+(x3-2)2 

/3(x)  =  xi  +x2  +  x3  -  1 

/4(x)  =  Xi  +  x2  -  x3  +  1 

/5(ar)  =  2x\  +  6x22  +  2(5x3  -  Xl  +  l)2 

/6(x)  =x\-  9x3 

Test  problem  8  (Madsen,  see  Overton  and  Murray,  1981) 

f\(x)  =  x\  +  xl  4-  xix2 
f2(x)  =  sin(xi) 
f3{x)  =  cos{x2) 


Test  problem  9  (Osborne  1,  1972) 

fi(x)  =  m  -  (zi  +  x2exp(-tlx4)  +  x3exp(-tiXs)) 

where  i  =  1,2,  •  •  •   ,33,  tz  =  10(i  —  1),  and 


i 

Vz 

i 

Vi 

i 

y% 

1 

0.844 

12 

0.718 

23 

0.478 

2 

0.908 

13 

0.685 

24 

0.467 

3 

0.932 

14 

0.658 

25 

0.457 

4 

0.936 

15 

0.628 

26 

0.448 

5 

0.925 

16 

0.603 

27 

0.438 

6 

0.908 

17 

0.580 

28 

0.431 

7 

0.881 

18 

0.558 

29 

0.424 

8 

0.850 

19 

0.538 

30 

0.420 

9 

0.818 

20 

0.522 

31 

0.414 

10 

0.784 

21 

0.506 

32 

0.411 

11 

0.751 

22 

0.490 

33 

0.406 

Test  problem  10  (Osborne  2) 


fl(x)  =yl  -  (xiexp(-ttx5)  +  x2exp{-(ti  -  x9)2x6) 

-\-x3exp(-{tt  -  x10)2a-7)  +  x4exp(-{ti  -  xu)2x$)) 


where  i  =  1,2,  •••   ,65,  tz  =  (i  —  1)/10,  and 


i 

Vi 

i 

Vi 

i 

Vr 

1 

1.366 

23 

0.694 

45 

0.672 

2 

1.191 

24 

0.644 

46 

0.708 

3 

1.112 

25 

0.624 

47 

0.633 

4 

1.013 

26 

0.661 

48 

0.668 

5 

0.991 

27 

0.612 

49 

0.645 

6 

0.885 

28 

0.558 

50 

0.632 

7 

0.831 

29 

0.533 

51 

0.591 

8 

0.847 

30 

0.495 

52 

0.559 

9 

0.786 

31 

0.500 

53 

0.597 

10 

0.725 

32 

0.423 

54 

0.625 

11 

0.746 

33 

0.395 

55 

0.739 

12 

0.679 

34 

0.375 

56 

0.710 

13 

0.608 

35 

0.372 

57 

0.729 

14 

0.655 

36 

0.391 

58 

0.720 

15 

0.616 

37 

0.396 

59 

0.636 

16 

0.606 

38 

0.405 

60 

0.581 

17 

0.602 

39 

0.428 

61 

0.428 

18 

0.626 

40 

0.429 

62 

0.292 

19 

0.651 

41 

0.523 

63 

0.162 

20 

0.724 

42 

0.562 

64 

0.098 

21 

0.649 

43 

0.607 

65 

0.054 

22 

0.649 

44 

0.653 

Test  problem  11  (Powell,  1962) 


/l(x)  =Xi  +  10X2 

/2(x)  =  51/2(.r3-^4) 
f3(x)  =  (x2-2x3)2 
f4(x)  =  l01/2(x1-x4)2 

Test  problem  12  (Rosenbrock,  1960) 

Mx)  =  I0(x2  -  x2,) 
f2(x)  =  1-  x1 

Test  problem  13  (Watson,  see  Kowalik  and  Osborne,  1968) 

fi(x)  =  Y^(j  -  l)xjt{-2  -  f^xj*}"1)     -1, 
where  i  =  1,  •  •  • ,  29,  U  =  i/29,  f3o(x)  =  x{  and  /3i(x)  =  x2  -  arf  -  1. 

Test  problem  14  (Wood,  see  Colville,  1968) 

fl(x)  =  10(x2-x2) 

f2(x)  =  1  -  Xi 

h(x)  =  901/2(x,-x2) 

/4(x)  =  1  -  x3 

f5(x)  =  l^2(x2+x4-2) 

f6(x)  =  10-1/2(x2-x4) 


Appendix  B 

"nlrq"    <- 

function(x,    y,    model,     t,    theta,    k    =    2,    eps    =    le-06,    big    =    le+20,    eta    =    0.97) 

{ 

#  This  is  a  function  to  compute  nonlinear  11  estimate. 

#  Input 

#  model  -  user-provided  function  which  returns  components 

#  f=(f_i  (x_i  ,  t) 

#  J=(grad  f_i  ) 

#  t      -  vector  of  initial  values  of  the  unknown  parameters 

#  theta  -  desired  quantile 

#  k      -  number  of  Meketon's  iterations  used  to  calculate  the 

#  dual  step 

#  eps    -  small  positive  number 

#  big    -  big  positive  number 

#  eta    -  0.97 

#  Output 

#  t      -  vector  of  estimated  parameters 

#  f        function  value  at  minimum 
# 

n  <-  length(y) 

zero  <-  rep ( 0 ,  n) 

w  <-  zero 

d  <-  rep(l,  n) 

m  <-  model (x,  y,  t,  theta) 

snew  <-  sum(abs (m$f ) ) 

sold  <-  big 

while (sold  -  snew  >  eps)  { 

z  <-  mekrq(m$J,  m$f,  w,  theta,  k,  int  =  F,  eps,  big,  eta) 

step  <-  z$coef 

#  Calculate  an  optimal  step  length  lambda 

lambda  <-  step . length ( t , step) 

t  <-  t  +  lambda  *  step 

m  <-  model (x,  y,  t,  theta) 

sold  <-  snew 

snew  <-  sum(abs (m$f ) ) 

w  <-  z$w 

w  <-  lsfit(m$J,  w,  int  =  F) $resid 

if (max (abs (w) )  >=  1) 

w  <-  w/ (max (abs (w) )  +  eps) 
} 

f  <-  snew 
return(t,  f) 
} 

"mekrq"  <- 

function(x,  y,  w,  theta,  kmax  =  1000,  int  =  T,  eps,  big,  eta) 

{ 

#  Compute  linear  regression  quantile  estimate  (Meketon,  1985). 

#  However,  note  that  W  is  not  initialized  and  the  maximum  number  of 

#  iteration  is  given  by  kmax. 
# 

if(int  ==  T)  x  <-  cbindd,  x) 

sr  <-  big 

k  <-  1 

while(k  <=  kmax  &  sr  -  crossprod(y,  w)  >  eps)  ( 

d  <-  pmin (theta  -  w,  1  -  theta  +  w) 

z  <-  lsfit(x,  y,  d~2,  int  =  F) 


sr  <-  sum(abs (z$resid) ) 

k  <-  k  +  1 

s  <-  z$resid  *  d~2 

alpha  <-  max(pmax(s/(theta  -  w) ,   -  s/(l  -  theta  +  w) 

w  <-  w  +  (eta/alpha)  *  s 

} 

coef  <-  z$coef 

returnfcoef ,  w) 


/* 

*  * 
** 
** 
** 

*  * 
** 
** 

*  * 
** 
** 
*• 
*• 
** 
** 
** 
** 
** 
** 
** 
** 
** 

* 
* 
* 
* 
* 

*/ 


/* 

* 

*/ 


This  procedure  is  a  function  to  compute  nonlinear  regression 
quantile  estimate. 

PROC  NLRQ 
FORMAT 

{  t,f  }  =  nlrq(  x,y, Smodel, start , theta, k, eps, big, eta  ) 
INPUT 

model  -  user-provided  function  which  returns  components 
f  =  (f_i  (x_i,t)) 
J  =  (grad  f_i) 
start  -  vector  of  initial  values  of  the  unknown  parameters 
theta  -  desired  quantile 

k  -  number  of  Meketon's  iterations  used  to  calcuate 
the  dual  step 
eps  -  small  positive  number  (le-06) 
big  -  big  positive  number  (le+20) 
eta  -  0.97 
OUTPUT 

t  -  vector  of  estimated  parameters 
f  -  function  value  at  minimum 


The  following  procedure  is  a  function  to  compute  linear 

regression  quantile  estimate  (Meketon,  1985).  However,  note 
that  W  is  not  initialized  and  the  maximum  number  of 
iterations  is  given  by  kmax. 

proc ( 2 )  =  mekrq(x,y,w, theta, kmax, inter, eps, big, eta) ; 
local  sr,k,wy,wx,d,t,r, s, alpha; 
if  inter  eq  1; 

x  =  ones (rows (x) , 1 ) -x; 
endif ; 
sr  =  big; 
k  =  1; 
do  while  k  <=  kmax  and  sr  -  y'*w  >  eps; 

d  =  minc( (theta-w) ' | ( 1-theta+w) ' ) ; 

d  =  vec(d) ; 

wx  =  x.*d; 

wy  =  y.*d; 

{  t  >  =  olsqr (wy,wx) ; 

r  =  y  -  x*t; 

sr  =  sumc (abs (r ) ) ; 

k  =  k  +  1; 

s  =  r.*dA2; 

alpha  =  maxc (maxc (( s. / (theta-w) )' | 
(-s./( 1-theta+w) ) ') ) ; 

w  =  w  +  eta/alpha*s; 
endo; 

retp( t,w) ; 
endp; 


proc  (  2  )  =  nlrq(x,y,  Smodel,  t,  theta,kmax,eps,big,eta)  ; 

local  n,d,w, f , J, snew, sold, step, lambda, t,  tl ,  p, inter, 

model :proc; 

n  =  rows (y ) ; 

w  =  zeros (n, 1 ) ; 

d  =  ones(n, 1 ) ; 

{  f,J  >  =  model (x,y,t) ; 

snew  =  sumc (abs (f ) ) ; 

sold  =  big; 

inter  =0; 

do  while  sold  -  snew  >  eps; 

{  step,w  }  =  mekrq( J, f ,w, theta, kmax, inter , eps, 
big, eta) ; 

/*  Calculate  an  optimal  step  length  lambda  */ 

{  lambda  }  =  stepl ( J, f ,t, step, &model,eps ) ; 

t  =  t  +  lambda*step; 

{  f,J  }  =  model (x,y,t) ; 

sold  =  snew; 

snew  =  sumc (abs ( f) ) ; 

{  tl,w,p  }  =  olsqr2(w,J); 

if  maxc(abs(w))  >=  1; 

w  =  w/ (maxc (abs (w) ) +eps) ; 

endif ; 
endo; 

retp( t, snew) ; 
endp; 
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